#include "Polynomial.h"

template<class Real>
Polynomial<Real>::Polynomial(int mdegree)
{
	if ( mdegree > 0 )
	{
		degree = mdegree;
		coeff = new Real[degree+1];
	}
	else
	{
		// default creation
		degree = -1;
		coeff = NULL;
	}
}

template <class Real>
Polynomial<Real>::Polynomial(const Polynomial<Real> &poly)
{
	degree = poly.degree;
	coeff = new Real[degree+1];
	for ( int i = 0; i < degree; i++ ) 
		coeff[i] = poly.coeff[i];
}

template <class Real>
Polynomial<Real>::~Polynomial()
{
	delete[] coeff;
}

template <class Real>
void Polynomial<Real>::set_degree(int degree)
{
	degree = mdegree;
	delete[] coeff;

	if ( degree >= 0 )
		coeff =  new Real[degree+1];
	else
		coeff = NULL;
}

template <class Real>
int Polynomial<Real>::get_degree() const
{
	return degree;
}

template <class Real>
Polynomial<Real>::operator const Real *() const
{
	return coeff;
}

template <class Real>
Polynomial<Real>::operator Real * ()
{
	return coeff;
}

template <class Real>
Real Polynomial<Real>::operator[] ( int i ) const
{
	assert( 0 <= i && i <= degree );
	return coeff[i];
}

template <class Real>
Polynomial<Real>& Polynomial<Real>::operator= ( const Polynomial& poly )
{
	delete[] coeff;
	degree = poly.degree;

	if ( degree >= 0 )
	{
		coeff = new Real[degree+1];
		for ( int i = 0; i <= degree; i++ )
			coeff[i] = poly.coeff[i];
	}

	return *this;
}


template <class Real>
Real Polynomial<Real>::operator () ( Real t ) const
{
	assert( degree >= 0 );

	Real result = coeff[degree];
	for ( int i = degree - 1; i >= 0; i-- )
	{
		result *= t;
		result += coeff[i];
	}
	return result;
}

template <class Real>
Polynomial<Real> Polynomial<Real>::operator+ ( const Polynomial& poly ) const
{
	assert ( degree >= 0 && poly.degree >= 0 );

	Polynomial<Real> sum;

	int i;

	if ( degree > poly.degree )
	{
		sum.set_degree( degree );
		for ( i = 0; i <= poly.degree; i++ )
			sum.coeff[i] = coeff[i] + poly.coeff[i];
		for ( i = poly.degree + 1; i <= degree; i++ )
			sum.coeff[i] = poly.coeff[i];
	}
	else
	{
		sum.set_degree( poly.degree );
		for ( i = 0; i <= degree; i++ )
			sum.coeff[i] = coeff[i] + poly.coeff[i];
		for ( i = degree + 1; i <= poly.degree; i++ )
			sum.coeff[i] = poly.coeff[i];
	}

	return sum;
}

template <class Real>
Polynomial<Real> Polynomial<Real>::operator- ( const Polynomial& poly ) const
{
	assert ( degree >= 0 && poly.degree >= 0 );

	Polynomial diff;
	int i;

	if ( degree > poly.degree )
	{
		diff.set_degree( degree );
		for ( i = 0 ; i <= poly.degree; i++ )
			diff.coeff[i] = coeff[i] - poly.coeff[i];
		for ( i = poly.degree+1; i <= degree; i++ )
			diff.coeff[i] = coeff[i];
	}
	else
	{
		diff.set_degree( poly.degree );
		for ( i = 0; i <= degree; i++ )
			diff.coeff[i] = coeff[i] - poly.coeff[i];
		for ( i = degree + 1; i <= poly.degree; i++ )
			diff.coeff[i] = -poly.coeff[i];
	}

	return diff;
}


template <class Real>
Polynomial<Real> Polynomial<Real>::operator* (const Polynomial<Real> &poly) const
{
	assert( degree >= 0 && poly.degree >= 0 );

	Polynomial<Real> prod( degree + poly.degree );

	memset( prod.coeff,0,(prod.degree+1)*sizeof(Real) );
	for ( int i0 = 0; i0 <= degree; i0++ )
	{
		for ( int i1 = 0; i1 <= poly.degree; i1++ )
		{ 
			int i2 = i0 + i1;
			prod.coeff[i2] += coeff[i0]*poly.coeff[i1];
		}
	}

	return prod;
}

template <class Real>
Polynomial<Real> Polynomial<Real>::operator+ (Real scalar) const
{
	assert( degree >= 0 );
	Polynomial sum( degree );
	sum.coeff[0] += scalar;
	return sum;
}

template <class Real>
Polynomial<Real> Polynomial<Real>::operator -(Real scalar) const
{
	assert( degree >= 0 );
	Polynomial diff( degree );
	diff.coeff[0] -= scalar;
	return diff;
}

template <class Real>
Polynomial<Real> Polynomial<Real>::operator /(Real scalar) const
{
	assert( degree >= 0 );
	Polynomial prod( degree );

	int i;
	if ( scalar != (Real) 0.0 )
	{
		Real inv_scalar = ( (Real) 1.0 )/scalar;
		for ( i = 0; i <= degree; i++ )
			prod.coeff[i] = inv_scalar * coeff[i];
	}
	else
	{
		for ( i = 0; i <= degree; i++ )
			prod.coeff[i] = Math<Real>::INFINITY;
	}

	return prod;
}

template <class Real>
Polynomial<Real> Polynomial<Real>::operator -()
{
	assert( degree >= 0 );
	Polynomial neg( degree );
	for ( int i = 0; i <= degree; i++ )
		neg.coeff[i] = -coeff[i];

	return neg;
}

template <class Real>
Polynomial<Real>& Polynomial<Real>::operator+= (const Polynomial<Real>& poly)
{
	assert( degree >= 0 );
	*this = *this + poly;
	return *this;
}

template <class Real>
Polynomial<Real>& Polynomial<Real>::operator-= (const Polynomial<Real>& poly)
{
	assert( degree >= 0 );
	*this = (*this) - poly;
	return *this;
}

template <class Real>
Polynomial<Real>& Polynomial<Real>::operator*= (const Polynomial<Real>& poly)
{
	assert( degree >= 0 );
	*this = (*this) * poly;
	return *this;
}

template <class Real>
Polynomial<Real>& Polynomial<Real>::operator+= (Real scalar)
{
	assert( degree >= 0 );
	coeff[0] += scalar;
	return *this;
}

template <class Real>
Polynomial<Real>& Polynomial<Real>::operator-= (Real scalar)
{
	assert( degree >= 0 );
	coeff[0] -= scalar;
	return *this;
}

template <class Real>
Polynomial<Real>& Polynomial<Real>::operator*= (Real scalar)
{
	assert( degree >= 0 );
	*this = (*this) * scalar;
	return *this;
}

template <class Real>
Polynomial<Real>& Polynomial<Real>::operator/= (Real scalar)
{
	assert( degree >= 0 );
	*this =(*this)/scalar;
	return *this;
}

template <class Real>
Polynomial<Real> Polynomial<Real>::derivative() const
{
	if ( degree > 0 )
	{
		Polynomial deriv( degree - 1 );
		for ( int i0 = 0; i1 = 1; i0 < degree; i0++; i1++ )
			deriv.coeff[i0] = i1 * coeff[i1];
		return deriv;
	}
	else if ( degree == 1 )
	{
		Polynomial mconst(0);
		mconst.coeff[0] = (Real) 0.0;
		return mconst;
	}
	return Polynomial();
}

template <class Real>
void Polynomial<Real>::compress( Real epsilon )
{
	int i;
	for ( i = degree; i >= 0; i-- )
	{
		if ( Math<Real>::abs( coeff[i] ) <= epsilon )
			degree--;
		else
			break;
	}

	if ( degree >= 0)
	{
		Real inv_leading = ((Real) 1.0)/coeff[degree];
		coeff[degree] = (Real)1.0;
		for ( i = 0; i < degree; i++ )
			coeff *= inv_leading;
	}
}

template <class Real>
void Polynomial<Real>::divide(const Polynomial<Real> &div, Polynomial<Real> &quot, Polynomial<Real> &rem, Real epsilon) const
{
	int Qdegree = degree - div.degree;
	if ( Qdegree >= 0 )
	{
		quot.set_degree( Qdegree );

		// temporary storage for the remainder
		Polynomial tmp = *this;

		// do the division ( Euclidean algorithm )
		Real inv = ((Real)1.0)/div[div.degree];
		for ( int iq = Qdegree; iq >= 0; iq-- )
		{
			int ir = div.degree + iq;
			quot[iq] = inv*tmp[ir];
			for ( ir--; ir >= iq; ir-- )
				tmp[ir] -= quot[iq]*div[ir-iq];
		}

		// calculate the correct degree for the remainder
		int rem_deg = div.degree - 1;
		while ( rem_deg > 0 && Math<Real>::abs( tmp[rem_deg] ) < epsilon )
			rem_deg--;

		if ( rem_deg == 0 && Math<Real>::abs( tmp[0] ) < epsilon )
			tmp[0] = ( Real ) 0.0;

		rem.set_degree(0);
		memcpy( rem.coeff, tmp.coeff, (rem_deg+1)*sizeof(Real) );
	}
	else
	{
		quot.set_degree(0);
		quot[0] = (Real)0.0;
		rem = *this;
	}
}

